function [P_music, f] = music(x, fs, num_sources, M, NFFT)
% MUSIC线谱估计函数（适用于复正弦信号）
% 输入参数：
%   x           : 输入信号（一维复向量）
%   fs          : 采样频率 (Hz)
%   num_sources : 信号源数目（频率成分数）
%   NFFT        : 频率点数（默认1024）
% 输出参数：
%   P_music     : MUSIC空间谱估计值（dB）
%   f           : 对应的频率向量 (Hz)

if nargin < 4
    NFFT = 1024; % 默认频率点数
end

% x = x(:); % 确保输入为列向量
N = length(x);

% 1. 计算自相关矩阵
% M = num_sources + 1; % 自相关矩阵阶数（需大于信号源数）
X = zeros(M, N-M+1);
for i = 1:N-M+1
    X(:, i) = x(i:i+M-1); % 构造数据矩阵（Hankel结构）
end
R = (X * X') / (N-M+1); % 样本协方差矩阵

% 2. 特征分解
[U, D] = eig(R); % U: 特征向量，D: 特征值
[~, idx] = sort(diag(D), 'descend'); % 特征值降序排序
U = U(:, idx);

% 提取噪声子空间
Un = U(:, num_sources+1:end); % 噪声子空间（小特征值对应的特征向量）

% 3. 计算频率网格
f = linspace(-fs/2, fs/2-1, NFFT); % 复信号需考虑负频率
P_music = zeros(NFFT, 1);

% 4. 计算MUSIC谱
for k = 1:NFFT
    a = transpose(exp(-1j * 2*pi * f(k)/fs * (0:M-1))); % 导向矢量（复指数）
    P_music(k) = 1 / (a' * (Un * Un') * a);   % 空间谱：1/(a^H Un Un^H a)
end

% 仅保留正频率部分（可选）
% P_music = fftshift(P_music); % 将负频率移到左侧
% f = fftshift(f);
% idx_pos = f >= 0;
% f = f(idx_pos);
% P_music = P_music(idx_pos);
end